CD ROM Paradise Collection 4/CD ROM Paradise Collection 4 1995 Nov.iso/graphics/3dvect37.zip/MATH.INC

Text File | 1994-06-22 | 83KB | 3,124 lines

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Filename     : Math.inc
; Included from: 3D1.ASM, 3D2.ASM, 3D3.ASM
; Description  : General math functions.
;
; Written by: John McCarthy
;             1316 Redwood Lane
;             Pickering, Ontario.
;             Canada, Earth, Milky Way (for those out-of-towners)
;             L1X 1C5
;
; Internet/Usenet:  [e-mail]
;         Fidonet:  Brian McCarthy 1:229/15
;   RIME/Relaynet: ->CRS
;
; Home phone, [phone], don't call at 2 am eh!
;
; Send me your protected mode source code!
; Send me your Objects!
; But most of all, Send me a postcard!!!!
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           public rotate        ; rotate using vmatrix
           public make3d        ; calculate 3d ?actual*?/z ( both x and y)
           public make3dx       ; xactual*x/z
           public make3dy       ; yactual*y/z
           public erotate       ; 32 bit rotate using ematrix
           public zsolve        ; solve single equation variable
           public ysolve
           public xsolve
           public cosign
           public sign
           public arctan
           public compound      ; generate rotation matrix (includes camera)
           public setsincose    ; set camera matrix
           public temp_matrix   ; set user defined/temporary matrix
           public temp_rotate   ; rotate point by temp matrix
           public matrix_multiply ; multiply tmatrix by vmatrix

           public set_precal7
           public set_precal147
           public frotate

           public fzsolve
           public fxsolve
           public fysolve

           public precal1
           public precal4
           public precal7

           public sqrt          ; eax=sqr(eax), thanks to TRAN!
           public _squareroot32 ; eax=sqr(eax)
           public sqrax2bx2     ; ax = sqr(ax^2+bx^2)

           public pre_cal_lambert   ; scan object si and calculate surface normals
           public calc_normal       ; guess...from 3 points, returns vector ebx,ecx,ebp
           public calc_d            ; calculate D from plane equation
           public lambert           ; calculate surface normal rotation matrix for object si
           public set_up_all_lambert; scans objects from si to di and calls pre_cal_lambert
           public lrotate           ; given normal for surface, figures out intensity

           public lx1               ; points to load up before calling calc_normal
           public ly1
           public lz1
           public lx2
           public ly2
           public lz2
           public lx3
           public ly3
           public lz3

           align 16

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Rotate - 32 bit rotate point using vmatrix
; In:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;     vmatrix - 32 bit rotation matrix - set up by "compound" routine
; Out:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;
; Notes:
;
; All rotations (erotate,rotate,temp_rotate) are 32 bit.
; Frotate uses rotation along a plane and uses ematrix with precal147
;
; point rotation
; ebx = x   ecx = y   ebp = z    32 bit rotation!
; clobbers edx,esi,eax
;
; remember , matrix offsets are:
;
;  0 1 2     multiply those by 4 for the word address of the matrix
;  3 4 5
;  6 7 8
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

rotate:
           mov eax,vmatrix+8    ; solve x = bx(0)+cx(1)+bp(2)
           imul ebp
           shrd eax,edx,14
           mov edi,eax
           mov eax,vmatrix+4
           imul ecx
           shrd eax,edx,14
           add edi,eax
           mov eax,vmatrix+0
           imul ebx
           shrd eax,edx,14
           add edi,eax   ; di = new x

           mov eax,vmatrix+20   ; solve y = bx(3)+cx(4)+bp(5)
           imul ebp
           shrd eax,edx,14
           mov esi,eax
           mov eax,vmatrix+16
           imul ecx
           shrd eax,edx,14
           add esi,eax
           mov eax,vmatrix+12
           imul ebx
           shrd eax,edx,14
           add esi,eax   ; si = new y

           mov eax,vmatrix+32   ; solve z = bx(6)+cx(7)+bp(8)
           imul ebp
           shrd eax,edx,14
           mov ebp,eax
           mov eax,vmatrix+28
           imul ecx
           shrd eax,edx,14
           add ebp,eax
           mov eax,vmatrix+24
           imul ebx
           shrd eax,edx,14
           add ebp,eax   ; bp = new z

           mov ecx,esi
           mov ebx,edi

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Make3d - scale 3d point into 2d point
; In:
;    EBX - x point
;    ECX - y point
;    EBP - z point
; Out:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;
; Notes:
;
; fast ratios found in macros.inc since
; multiplication has been substituted with fast lea
;
; trashes eax,edx,edi
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

make3d:                            ; bp must always be positive
           cmul eax,ebx,ratiox     ; use fast constant multiply

           idiv ebp
           mov ebx,eax

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Make3dy - scale 3d point into 2d point on x axis only
; In:
;    ECX - y point
;    EBP - z point
; Out:
;    ECX - y point
;    EBP - z point
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
make3dy:
           cmul eax,ecx,ratioy

           idiv ebp
           mov ecx,eax

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Make3dx - scale 3d point into 2d point on y axis only
; In:
;    EDI - x point
;    ESI - z point
; Out:
;    EDI - x point
;    ESI - z point
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

make3dx:                           ; bp must always be positive
           cmul eax,edi,ratiox

           idiv esi
           mov edi,eax

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Checkfront: checks if a side is visible. (counter-clockwise)
;
; In:
;   (EDI,EBP) - xy of point 1
;   (ESI,ECX) - xy of point 2
;   (EDX,EBX) - xy of point 3
; Out:
;   ECX < 0 if side counter-clockwise
;
; Notes: routine courtesy of "RAZOR"
; eg:
;          call checkfront
;          cmp ecx,0
;          jng dontdraw
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░


           align 16

checkfront:
           cmp edi,esi
           jng s cfc2
           mov eax,edi
           mov edi,esi
           mov esi,edx
           mov edx,eax
           mov eax,ebp
           mov ebp,ecx
           mov ecx,ebx
           mov ebx,eax
cfc:
           cmp edi,esi
           jng s cfc2
           mov eax,edi
           mov edi,esi
           mov esi,edx
           mov edx,eax
           mov eax,ebp
           mov ebp,ecx
           mov ecx,ebx
           mov ebx,eax
cfc2:
           mov eax,edx               ; ax = x3
           sub eax,edi               ; ax = x3 - x1
           sub ecx,ebp               ; cx = y2 - y1
           imul ecx                  ; ax = (x3-x1)*(y2-y1)
           mov ecx,eax               ; save it...
           mov eax,esi               ; ax = x2
           sub eax,edi               ; ax = x2 - x1
           sub ebx,ebp               ; bx = y3 - y1
           imul ebx                  ; ax = (x2-x1)*(y3-y1)
           sub ecx,eax               ; cx = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1)
           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; ERotate - 32 bit rotate point using ematrix
; In:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;     ematrix - 32 bit rotation matrix - set up by "setsincose" routine
; Out:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;
; Notes:
;
; point rotation for eye - solves all x,y,z parameters
; camera rotation is 32 bit and uses ematrix
;
; remember , matrix offsets are:
;
;  0 1 2     multiply those by 4 for the doubleword address of the matrix
;  3 4 5
;  6 7 8
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           align 16

erotate:
           mov eax,ematrix+8
           imul ebp
           shrd eax,edx,14
           mov edi,eax
           if usez eq yes
           mov eax,ematrix+4
           imul ecx
           shrd eax,edx,14
           add edi,eax
           endif
           mov eax,ematrix+0
           imul ebx
           shrd eax,edx,14
           add edi,eax   ; di = new x

           mov eax,ematrix+20
           imul ebp
           shrd eax,edx,14
           mov esi,eax
           mov eax,ematrix+16
           imul ecx
           shrd eax,edx,14
           add esi,eax
           mov eax,ematrix+12
           imul ebx
           shrd eax,edx,14
           add esi,eax   ; si = new y

           mov eax,ematrix+32
           imul ebp
           shrd eax,edx,14
           mov ebp,eax
           mov eax,ematrix+28
           imul ecx
           shrd eax,edx,14
           add ebp,eax
           mov eax,ematrix+24
           imul ebx
           shrd eax,edx,14
           add ebp,eax   ; bp = new z

           mov ecx,esi
           mov ebx,edi

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Zsolve - 32 bit rotate point using ematrix - solve one variable only
; In:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;    ematrix - 32 bit rotation matrix - set up by "setsincose" routine
;
; Out:
;    EBX - x point (same as entry)
;    ECX - y point (same as entry)
;    EBP - z point (same as entry)
;    ESI - new z point/location
;
; Notes:
;
; solve z from ematrix - same as above erotate but only solves z for fast
; test of where object is - result is in esi
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           align 16

zsolve:
           mov eax,ematrix+32
           imul ebp
           shrd eax,edx,14
           mov esi,eax
           mov eax,ematrix+28
           imul ecx
           shrd eax,edx,14
           add esi,eax
           mov eax,ematrix+24
           imul ebx
           shrd eax,edx,14
           add esi,eax   ; si = new z
           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Xsolve - 32 bit rotate point using ematrix - solve one variable only
; In:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;    ematrix - 32 bit rotation matrix - set up by "setsincose" routine
;
; Out:
;    EBX - x point (same as entry)
;    ECX - y point (same as entry)
;    EBP - z point (same as entry)
;    EDI - new x point/location
;
; Notes:
; if object z test from above routine is positive, this routine will solve
; the rest of the rotation matrix.  this is so we don't waste time solving
; for x and y locations if the object is behind the camera anyway.
; saves imuls
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           align 16
xsolve:
           mov eax,ematrix+8
           imul ebp
           shrd eax,edx,14
           mov edi,eax
           if usez eq yes
           mov eax,ematrix+4
           imul ecx
           shrd eax,edx,14
           add edi,eax
           endif
           mov eax,ematrix+0
           imul ebx
           shrd eax,edx,14
           add edi,eax   ; di = new x
           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Ysolve - 32 bit rotate point using ematrix - solve one variable only
; In:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;    ESI - new z point
;    EDI - new x point
;    ematrix - 32 bit rotation matrix - set up by "setsincose" routine
;
; Out:
;    EBX - x new point from EDI
;    ECX - y new point
;    EBP - z new point from ESI
;
; Notes:
;
; solve y from ematrix - same as above xsolve but solves y for fast
; test of where object is.  final variables are then cleaned up to
; immitate the erotate function in parts.
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           align 16
ysolve:
           mov eax,ematrix+16
           imul ecx
           shrd eax,edx,14
           mov ecx,eax
           mov eax,ematrix+12
           imul ebx
           shrd eax,edx,14
           add ecx,eax
           mov eax,ematrix+20
           imul ebp
           shrd eax,edx,14
           add ecx,eax   ; cx = new y

           mov ebx,edi   ; final test, move into appropriate regs
           mov ebp,esi

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
;   Sign - 16 bit theta to 32bit sin(@)
; In:
;     AX - theta  0 - 65536 (0-360)
; Out:
;    EAX - sin (@)   (-4000h to 4000h)
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; CoSign - 16 bit theta to 32bit cos(@)
; In:
;     AX - theta  0 - 65536 (0-360)
; Out:
;    EAX - cos (@)   (-4000h to 4000h)
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Notes:
; calculate sin into eax, from ax, smashes bx
; after imul by sin, shr eax,14 to compensate for decimal factor!
;  eg:
;    mov eax,sin(@)
;    mov ebx,32bitnumber
;    imul ebx
;    shrd eax,edx,14
;    eax = ebx*sin(@)
;
;    mov ax,sin(@)
;    mov bx,16bitnumber
;    imul bx
;    shrd ax,dx,14
;    eax = bx*sin(@)
;
; eax is simply a sign extended ax
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           align 16

cosign:
           add ax,4000h
sign:
           shr ax,2
           cmp ax,2000h
           jge s q3o4         ; quadrant 3 or 4

           cmp ax,1000h
           jl s q0            ; quad 1

           mov ebx,1fffh
           sub bx,ax
           jmp s halfsign     ; quad 2
q0:
           movzx ebx,ax
           jmp s halfsign
q3o4:
           cmp ax,3000h
           jl s q3
           mov ebx,3fffh
           sub bx,ax
           call halfsign      ; quad 4
           neg eax
           ret
q3:
           and ax,0fffh
           movzx ebx,ax       ; quad 3
           call halfsign
           neg eax
           ret
halfsign:
           xor eax,eax
           mov ax,w sinus[ebx*2]
           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Arctan - 32 bit rise/run to 16bit arctan(rise/run)
; In:
;    EAX - Run
;    ECX - Rise
; Out:
;     AX - arctan(ECX/EAX)
;
; Notes:
; smashes cx,ax,dx,si
; arctan(ecx/0) is valid and tested for
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           align 16

arctan:
           cmp eax,0
           jl s qd2or3
           cmp ecx,0
           jge s halftax      ; quadrant 1
           neg ecx            ; quadrant 4, ax=-ax
           call halftan
           neg ax
           shl eax,2
           ret
qd2or3:
           neg eax
           cmp ecx,0
           jge s qd2
           neg ecx            ; quad 3, ax=ax+8192
           call halftan
           add ax,8192
           shl eax,2
           ret
qd2:
           call halftan
           neg ax
           add ax,8192
           shl eax,2
           ret
halftax:
           call halftan
           shl eax,2
           ret

           align 16

halftan:
           xor edx,edx

; cx=rise  positive
; ax=run   positive

           cmp eax,ecx
           jl s opptan        ; greater than 45 degrees, other side...

           xchg ecx,eax       ; ax<cx
           shld edx,eax,11    ; *2048 edx = high dword for divide
           shl eax,11         ; *2048
           div ecx
           movzx esi,ax
           mov ax,w negtan[esi*2] ; resulting angle (0-512 is 0-45) in ax
           ret

           align 16

opptan:
           shld edx,eax,11    ; *2048 edx = high dword for divide
           shl eax,11         ; *2048

           div ecx
           movzx esi,ax       ; ax remainder
           mov cx,w negtan[esi*2]
           mov eax,1000h
           sub ax,cx          ; resulting angle (2048-4096 is 45-90) in ax
           ret

           align 16

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Compound - generate object matrix, 12 imul's first
; In:
;    ESI - Object # to get angles from
;    vxs[esi*2] - object x angle (0-65536)
;    vys[esi*2] - object y angle (0-65536)
;    vzs[esi*2] - object z angle (0-65536)
; Out:
;    vmatrix - resulting rotation matrix including camera matrix
;    ESI = ESI
;
; Notes:
;              x                         y                      z
;
;x=  cz * cy - sx * sy * sz   - sz * cy - sx * sy * cz     - cx * sy
;
;y=         sz * cx                   cx * cz                - sx
;
;z=  cz * sy + sx * sz * cy   - sy * sz + sx * cy * cz       cx * cy
;
;then perform matrix multiply by negative x and z matricies
;
; -x matrix                             -z matrix
;     x       y       z                   x       y       z
;
;x    1       0       0                cz     sz       0
;
;y    0      cx       sx              -sz     cz       0
;
;z    0     -sx       cx                0      0       1
;
; notice original object matrix takes 12 imuls, camera modify takes 24, can
; you do this faster? (less imuls)
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

compound:
           push esi

           mov ax,vxs[esi*2]
           neg ax
           push eax
           call cosign
           mov vcosx,eax
           pop eax
           call sign
           mov vsinx,eax
           mov ebp,eax            ; bp = sx
           neg eax
           mov [vmatrix+20],eax

           mov ax,vzs[esi*2]
           neg ax
           push eax
           call cosign
           mov vcosz,eax
           mov edi,eax            ; di = cz
           pop eax
           call sign
           mov vsinz,eax
           mov edx,eax            ; dx = sz

           mov ax,vys[esi*2]
           neg ax
           add ax,eyeay
           push eax
           call cosign
           mov vcosy,eax
           pop eax
           call sign
           mov vsiny,eax          ; ax = sy

           mov ebx,edx            ; save sz

           mov ecx,eax            ; save sy

           imul ebx               ; bx = - sy * sz
           shr eax,14
           movsx ebx,ax
           neg ebx
           mov [vmatrix+28],ebx

           mov eax,ecx            ; si = cz * sy
           imul edi
           shr eax,14
           movsx esi,ax
           mov [vmatrix+24],esi

           mov eax,vcosy

           imul edi               ; di = cy * cz
           shr eax,14
           movsx edi,ax
           mov [vmatrix+0],edi

           mov eax,vsinz
           mov ecx,vcosy

           imul ecx               ; cx = - sz * cy
           shr eax,14
           movsx ecx,ax
           neg ecx
           mov [vmatrix+4],ecx

           mov eax,ebp
           imul esi
           shr eax,14
           movsx esi,ax
           neg esi
           add [vmatrix+4],esi

           mov eax,ebp
           imul edi
           shr eax,14
           movsx edi,ax
           add [vmatrix+28],edi

           mov eax,ebp
           imul ebx
           shr eax,14
           movsx ebx,ax
           add [vmatrix+0],ebx

           mov eax,ebp
           imul ecx
           shr eax,14
           movsx ecx,ax
           neg ecx
           add [vmatrix+24],ecx

           mov esi,vcosx

           mov eax,vcosy
           imul esi                   ; cx * cy
           shr eax,14
           movsx eax,ax
           mov [vmatrix+32],eax

           mov eax,vsiny
           imul esi                   ;-cx * sy
           shr eax,14
           movsx eax,ax
           neg eax
           mov [vmatrix+8],eax

           mov eax,vsinz
           imul esi                   ; cx * sz
           shr eax,14
           movsx eax,ax
           mov [vmatrix+12],eax

           mov eax,vcosz
           imul esi                   ; cx * cz
           shr eax,14
           movsx eax,ax
           mov [vmatrix+16],eax

           mov edi,ecosx              ; now perform camera x rotation,12 imuls
           mov esi,esinx
           mov ebp,esi
           neg ebp

           mov eax,[vmatrix+12]
           imul edi
           shr eax,14
           movsx ecx,ax

           mov eax,[vmatrix+24]
           imul esi
           shr eax,14
           movsx eax,ax

           add ecx,eax                ; ecx = new vmatrix+12

           mov eax,[vmatrix+12]
           imul ebp
           shr eax,14
           movsx ebx,ax

           mov eax,[vmatrix+24]
           imul edi
           shr eax,14
           movsx eax,ax

           add ebx,eax                ; ebx = new vmatrix+24

           mov [vmatrix+12],ecx
           mov [vmatrix+24],ebx

           mov eax,[vmatrix+16]
           imul edi
           shr eax,14
           movsx ecx,ax

           mov eax,[vmatrix+28]
           imul esi
           shr eax,14
           movsx eax,ax

           add ecx,eax                ; ecx = new vmatrix+16

           mov eax,[vmatrix+16]
           imul ebp
           shr eax,14
           movsx ebx,ax

           mov eax,[vmatrix+28]
           imul edi
           shr eax,14
           movsx eax,ax

           add ebx,eax                ; ebx = new vmatrix+28

           mov [vmatrix+16],ecx
           mov [vmatrix+28],ebx

           mov eax,[vmatrix+20]
           imul edi
           shr eax,14
           movsx ecx,ax

           mov eax,[vmatrix+32]
           imul esi
           shr eax,14
           movsx eax,ax

           add ecx,eax                ; ecx = new vmatrix+20

           mov eax,[vmatrix+20]
           imul ebp
           shr eax,14
           movsx ebx,ax

           mov eax,[vmatrix+32]
           imul edi
           shr eax,14
           movsx eax,ax

           add ebx,eax                ; ebx = new vmatrix+32

           mov [vmatrix+20],ecx
           mov [vmatrix+32],ebx

           if usez eq yes

           mov edi,ecosz              ; now perform camera z rotation,12 imuls
           mov esi,esinz
           mov ebp,esi
           neg esi

           mov eax,[vmatrix+0]
           imul edi
           shr eax,14
           movsx ecx,ax

           mov eax,[vmatrix+12]
           imul esi
           shr eax,14
           movsx eax,ax

           add ecx,eax

           mov eax,[vmatrix+0]
           imul ebp
           shr eax,14
           movsx ebx,ax

           mov eax,[vmatrix+12]
           imul edi
           shr eax,14
           movsx eax,ax

           add ebx,eax

           mov [vmatrix+0],ecx
           mov [vmatrix+12],ebx

           mov eax,[vmatrix+4]
           imul edi
           shr eax,14
           movsx ecx,ax

           mov eax,[vmatrix+16]
           imul esi
           shr eax,14
           movsx eax,ax

           add ecx,eax

           mov eax,[vmatrix+4]
           imul ebp
           shr eax,14
           movsx ebx,ax

           mov eax,[vmatrix+16]
           imul edi
           shr eax,14
           movsx eax,ax

           add ebx,eax

           mov [vmatrix+4],ecx
           mov [vmatrix+16],ebx

           mov eax,[vmatrix+8]
           imul edi
           shr eax,14
           movsx ecx,ax

           mov eax,[vmatrix+20]
           imul esi
           shr eax,14
           movsx eax,ax

           add ecx,eax

           mov eax,[vmatrix+8]
           imul ebp
           shr eax,14
           movsx ebx,ax

           mov eax,[vmatrix+20]
           imul edi
           shr eax,14
           movsx eax,ax

           add ebx,eax

           mov [vmatrix+8],ecx
           mov [vmatrix+20],ebx

           endif

           pop esi
           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Setsincose - generate rotation matrix for  y,x,z  camera rotation
;
; In:
;    eyeax - camera x angle (0-65536)
;    eyeay - camera y angle (0-65536)
;    eyeaz - camera z angle (0-65536)
; Out:
;    vmatrix - resulting rotation matrix including camera matrix
;
; Notes:
; called only once every frame.  completed in 12 multiplys
; matrix is also used for objects with no rotation (always angle 0,0,0)
;
; where is my postcard! see readme.doc for info.
;
;              x                    y                    z
;
; x=  cz * cy + sx * sy * sz     -cx * sz     - sy * cz + sx * cy * sz
;
; y=  sz * cy - sx * sy * cz      cx * cz     - sy * sz - sz * cy * cz
;
; z=         cx * sy                 sx                cx * cy
;
;
;  matrix offsets: (doublewords)
;
;     x  y  z
;
; x    0  4  8
; y   12 16 20
; z   24 28 32
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           align 16

setsincose:

           mov ax,eyeax
           call cosign
           mov ecosx,eax          ; ecosx and such are used by object rotation
           mov ax,eyeax           ; ematrix is used to find where object is
           call sign
           mov esinx,eax
           mov [ematrix+28],eax
           mov ebp,eax            ; bp = sx

           if usez eq yes
           mov ax,eyeaz
           call cosign
           mov ecosz,eax
           mov edi,eax            ; di = cz
           mov ax,eyeaz
           call sign
           mov esinz,eax
           mov edx,eax            ; dx = sz
           endif

           if usez eq no
           mov edi,4000h          ; di = cos 0
           mov ecosz,4000h
           xor edx,edx            ; dx = sin 0
           mov esinz,0
           endif

           mov ax,eyeay
           call cosign
           mov ecosy,eax
           mov ax,eyeay
           call sign
           mov esiny,eax          ; ax = sy

           mov ebx,edx            ; save sz

           mov ecx,eax            ; save sy

           imul bx                ; bx = sy * sz
           shrd ax,dx,14
           movsx ebx,ax
           neg ebx
           mov [ematrix+20],ebx
           neg ebx

           mov eax,ecx            ; si = - (cz * sy)
           imul di
           shrd ax,dx,14
           movsx esi,ax
           neg esi
           mov [ematrix+8],esi

           mov eax,ecosy

           imul di                ; di = cy * cz
           shrd ax,dx,14
           movsx edi,ax
           mov [ematrix+0],edi

           mov eax,esinz
           mov ecx,ecosy

           imul cx                ; cx = sz * cy
           shrd ax,dx,14
           movsx ecx,ax
           mov [ematrix+12],ecx

           mov eax,ebp
           imul si
           shrd ax,dx,14
           movsx esi,ax
           add [ematrix+12],esi

           mov eax,ebp
           imul di
           shrd ax,dx,14
           movsx edi,ax
           neg edi
           add [ematrix+20],edi

           mov eax,ebp
           imul bx
           shrd ax,dx,14
           movsx ebx,ax
           add [ematrix+0],ebx

           mov eax,ebp
           imul cx
           shrd ax,dx,14
           movsx ecx,ax
           add [ematrix+8],ecx

           mov esi,ecosx

           mov eax,ecosy
           imul si                    ; cx * cy
           shrd ax,dx,14
           movsx eax,ax
           mov [ematrix+32],eax

           mov eax,esiny
           imul si                    ; cx * sy
           shrd ax,dx,14
           movsx eax,ax
           mov [ematrix+24],eax

           mov eax,esinz
           imul si                    ;-cx * sz
           shrd ax,dx,14
           movsx eax,ax
           neg eax
           mov [ematrix+4],eax

           mov eax,ecosz
           imul si                    ; cx * cz
           shrd ax,dx,14
           movsx eax,ax
           mov [ematrix+16],eax

           neg esinx                  ; reverse angles for object rotation
           neg esiny

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Temp_Matrix: generate temp matrix, 12 imul's, from object esi
;
; In:
;    ESI - Object # to get angles from
;    vxs[esi*2] - object x angle (0-65536)
;    vys[esi*2] - object y angle (0-65536)
;    vzs[esi*2] - object z angle (0-65536)
; Out:
;    tmatrix - resulting rotation matrix (excluding camera matrix)
;    ESI = ESI
;
; Notes:
;              x                         y                      z
;
;x=  cz * cy - sx * sy * sz   - sz * cy - sx * sy * cz     - cx * sy
;
;y=         sz * cx                   cx * cz                - sx
;
;z=  cz * sy + sx * sz * cy   - sy * sz + sx * cy * cz       cx * cy
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

temp_matrix:
           push esi

           mov ax,vxs[esi*2]
           neg ax
           push eax
           call cosign
           mov vcosx,eax
           pop eax
           call sign
           mov vsinx,eax
           mov ebp,eax            ; bp = sx
           neg eax
           mov [tmatrix+20],eax

           mov ax,vzs[esi*2]
           neg ax
           push eax
           call cosign
           mov vcosz,eax
           mov edi,eax            ; di = cz
           pop eax
           call sign
           mov vsinz,eax
           mov edx,eax            ; dx = sz

           mov ax,vys[esi*2]
           neg ax
           push eax
           call cosign
           mov vcosy,eax
           pop eax
           call sign
           mov vsiny,eax          ; ax = sy

           mov ebx,edx            ; save sz

           mov ecx,eax            ; save sy

           imul ebx               ; bx = - sy * sz
           shr eax,14
           movsx ebx,ax
           neg ebx
           mov [tmatrix+28],ebx

           mov eax,ecx            ; si = cz * sy
           imul edi
           shr eax,14
           movsx esi,ax
           mov [tmatrix+24],esi

           mov eax,vcosy

           imul edi               ; di = cy * cz
           shr eax,14
           movsx edi,ax
           mov [tmatrix+0],edi

           mov eax,vsinz
           mov ecx,vcosy

           imul ecx               ; cx = - sz * cy
           shr eax,14
           movsx ecx,ax
           neg ecx
           mov [tmatrix+4],ecx

           mov eax,ebp
           imul esi
           shr eax,14
           movsx esi,ax
           neg esi
           add [tmatrix+4],esi

           mov eax,ebp
           imul edi
           shr eax,14
           movsx edi,ax
           add [tmatrix+28],edi

           mov eax,ebp
           imul ebx
           shr eax,14
           movsx ebx,ax
           add [tmatrix+0],ebx

           mov eax,ebp
           imul ecx
           shr eax,14
           movsx ecx,ax
           neg ecx
           add [tmatrix+24],ecx

           mov esi,vcosx

           mov eax,vcosy
           imul esi                   ; cx * cy
           shr eax,14
           movsx eax,ax
           mov [tmatrix+32],eax

           mov eax,vsiny
           imul esi                   ;-cx * sy
           shr eax,14
           movsx eax,ax
           neg eax
           mov [tmatrix+8],eax

           mov eax,vsinz
           imul esi                   ; cx * sz
           shr eax,14
           movsx eax,ax
           mov [tmatrix+12],eax

           mov eax,vcosz
           imul esi                   ; cx * cz
           shr eax,14
           movsx eax,ax
           mov [tmatrix+16],eax

           pop esi
           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Temp_Rotate - 32 bit rotate point using tmatrix
; In:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;     tmatrix - 32 bit rotation matrix - set up by "temp_matrix" routine
; Out:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;
; Notes:
;  Same as rotate and erotate
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

temp_rotate:
           mov eax,tmatrix+8    ; solve x = bx(0)+cx(1)+bp(2)
           imul ebp
           shrd eax,edx,14
           mov edi,eax
           mov eax,tmatrix+4
           imul ecx
           shrd eax,edx,14
           add edi,eax
           mov eax,tmatrix+0
           imul ebx
           shrd eax,edx,14
           add edi,eax   ; di = new x

           mov eax,tmatrix+20   ; solve y = bx(3)+cx(4)+bp(5)
           imul ebp
           shrd eax,edx,14
           mov esi,eax
           mov eax,tmatrix+16
           imul ecx
           shrd eax,edx,14
           add esi,eax
           mov eax,tmatrix+12
           imul ebx
           shrd eax,edx,14
           add esi,eax   ; si = new y

           mov eax,tmatrix+32   ; solve z = bx(6)+cx(7)+bp(8)
           imul ebp
           shrd eax,edx,14
           mov ebp,eax
           mov eax,tmatrix+28
           imul ecx
           shrd eax,edx,14
           add ebp,eax
           mov eax,tmatrix+24
           imul ebx
           shrd eax,edx,14
           add ebp,eax   ; bp = new z

           mov ecx,esi
           mov ebx,edi

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Matrix_Multiply: multiply tmatrix by vmatrix, [vmatrix]=[tmatrix][vmatrix]
;
; In:
;    vmatrix - rotation matrix
;    tmatrix - rotation matrix
; Out:
;    vmatrix - resulting rotation matrix
;
; Notes:
;
; [ tmatrix+ 0 tmatrix+ 2 tmatrix+ 4 ] [ vmatrix+ 0 vmatrix+ 2 vmatrix+ 4 ]
; [                                  ] [                                  ]
; [ tmatrix+ 6 tmatrix+ 8 tmatrix+10 ] [ vmatrix+ 6 vmatrix+ 8 vmatrix+10 ]
; [                                  ] [                                  ]
; [ tmatrix+12 tmatrix+14 tmatrix+16 ] [ vmatrix+12 vmatrix+14 vmatrix+16 ]
;
; Think of it this way, this routine will generate a resulting matrix as  if
; you rotated an  object by tmatrix, then rotated  the  object  by  vmatrix.
; Instead, call this routine then you will only have to rotate the object by
; Vmatrix!.
;
; Notice Tmatrix is done before Vmatrix!!  This is used for calculating the
; positions of arms on bodies, hands on arms, fingers on hands...etc...
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

matrix_multiply:

           mov ebx,[vmatrix+0]
           mov ecx,[vmatrix+4]
           mov ebp,[vmatrix+8]

           mov eax,[tmatrix+0]
           imul ebx
           shrd eax,edx,14
           mov esi,eax

           mov eax,[tmatrix+12]
           imul ecx
           shrd eax,edx,14
           add esi,eax

           mov eax,[tmatrix+24]
           imul ebp
           shrd eax,edx,14
           add esi,eax

           push esi             ; tmatrix+0

           mov eax,[tmatrix+4]
           imul ebx
           shrd eax,edx,14
           mov esi,eax

           mov eax,[tmatrix+16]
           imul ecx
           shrd eax,edx,14
           add esi,eax

           mov eax,[tmatrix+28]
           imul ebp
           shrd eax,edx,14
           add esi,eax

           push esi             ; tmatrix+4

           mov eax,[tmatrix+8]
           imul ebx
           shrd eax,edx,14
           mov esi,eax

           mov eax,[tmatrix+20]
           imul ecx
           shrd eax,edx,14
           add esi,eax

           mov eax,[tmatrix+32]
           imul ebp
           shrd eax,edx,14
           add esi,eax

           push esi             ; tmatrix+8

           mov ebx,[vmatrix+12]
           mov ecx,[vmatrix+16]
           mov ebp,[vmatrix+20]

           mov eax,[tmatrix+0]
           imul ebx
           shrd eax,edx,14
           mov esi,eax

           mov eax,[tmatrix+12]
           imul ecx
           shrd eax,edx,14
           add esi,eax

           mov eax,[tmatrix+24]
           imul ebp
           shrd eax,edx,14
           add esi,eax

           push esi             ; tmatrix+12

           mov eax,[tmatrix+4]
           imul ebx
           shrd eax,edx,14
           mov esi,eax

           mov eax,[tmatrix+16]
           imul ecx
           shrd eax,edx,14
           add esi,eax

           mov eax,[tmatrix+28]
           imul ebp
           shrd eax,edx,14
           add esi,eax

           push esi             ; tmatrix+16

           mov eax,[tmatrix+8]
           imul ebx
           shrd eax,edx,14
           mov esi,eax

           mov eax,[tmatrix+20]
           imul ecx
           shrd eax,edx,14
           add esi,eax

           mov eax,[tmatrix+32]
           imul ebp
           shrd eax,edx,14
           add esi,eax

           push esi             ; tmatrix+20

           mov ebx,[vmatrix+24]
           mov ecx,[vmatrix+28]
           mov ebp,[vmatrix+32]

           mov eax,[tmatrix+0]
           imul ebx
           shrd eax,edx,14
           mov esi,eax

           mov eax,[tmatrix+12]
           imul ecx
           shrd eax,edx,14
           add esi,eax

           mov eax,[tmatrix+24]
           imul ebp
           shrd eax,edx,14
           add esi,eax

           push esi             ; tmatrix+24

           mov eax,[tmatrix+4]
           imul ebx
           shrd eax,edx,14
           mov esi,eax

           mov eax,[tmatrix+16]
           imul ecx
           shrd eax,edx,14
           add esi,eax

           mov eax,[tmatrix+28]
           imul ebp
           shrd eax,edx,14
           add esi,eax

           push esi             ; tmatrix+28

           mov eax,[tmatrix+8]
           imul ebx
           shrd eax,edx,14
           mov esi,eax

           mov eax,[tmatrix+20]
           imul ecx
           shrd eax,edx,14
           add esi,eax

           mov eax,[tmatrix+32]
           imul ebp
           shrd eax,edx,14
           add esi,eax

;          push esi             ; tmatrix+32

;          pop esi
           mov [vmatrix+32],esi
           pop esi
           mov [vmatrix+28],esi
           pop esi
           mov [vmatrix+24],esi
           pop esi
           mov [vmatrix+20],esi
           pop esi
           mov [vmatrix+16],esi
           pop esi
           mov [vmatrix+12],esi
           pop esi
           mov [vmatrix+ 8],esi
           pop esi
           mov [vmatrix+ 4],esi
           pop esi
           mov [vmatrix+ 0],esi
           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Getroot:  ax = sqr(ax)
;
; In:
;    AX - number to take square root of (0 - 65535)
; Out:
;    AX - Square root
;
; Notes:
; I like TRANs method better - greater resolution and it doesn't require the
; table. Im not too sure whos method is faster.
;
; This routine is not used in the main animation loop.  It could be used for
; things like determining collision detection since collision does not require
; much accuracy and collision does require many many checks.
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
;           include squares.inc
;
;getroot:                      ; get square root of ax, where ax = 0-65535
;           cmp ax,0fe01h      ; since ax cannot be negative anyway!
;           jae sqr255         ; routine requires squares tables.
;           mov si,offset squares
;           mov cx,ax
;           inc cx
;           cld
;nextroot:
;           lodsw
;           cmp ax,cx
;           jbe  nextroot      ; jb is exact but jbe is better approximation
;           mov ax,si
;           sub ax,offset squares+3
;           sar ax,1
;           ret
;sqr255:
;           mov ax,255
;           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; _Squareroot32:  ax = sqr(ax)
; In:
;   EAX - number to take square root of (32bit)
; Out:
;   EAX = EBX = Square root
; Notes:
;  This came from Dave Stamp's VR386 - Hey, he's short, fat, and stupid looking
;  but I liked his routine anyway...
;
;  takes root of 32 bit number to 16 bit result
;  about 220 clocks worst case:
;  3 us on 486/66 and 10 us on 386/25
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

sqrt32  macro
        local skip
        shld edx,eax,2 ; get 2 bits of input to error
        shl eax,2
        add ebx,ebx    ; estimate*2
        mov ecx,ebx    ; temp = est*2
        add ecx,ecx
        cmp edx,ecx    ; error>2*est?
        jle s skip
        inc ebx        ; yes, update for new bit added
        inc ecx
        sub edx,ecx
skip:
        endm

        public _squareroot32

_squareroot32:
        test eax,0ffff0000h  ; can we cut it in half?
        jne s hasupper

        shl eax,16           ; yes, so prescale
        test eax,0ff000000h  ; half again?
        jne do16

        shl eax,8            ; yes, prescale
        jmp do8              ; do 8 loops

hasupper:
        test eax,0ff000000h  ; half again?
        jne do32
        shl eax,8
        jmp do24
do32:
        sqrt32
        sqrt32
        sqrt32
        sqrt32
do24:
        sqrt32
        sqrt32
        sqrt32
        sqrt32
do16:
        sqrt32
        sqrt32
        sqrt32
        sqrt32
do8:
        sqrt32
        sqrt32
        sqrt32
        sqrt32

        mov eax,ebx        ; result in eax and ebx
        ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Sqrt: Routine courtesy TRAN
;
; In:
;   EAX - number to take root of
; Out:
;   EAX - root
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
sqrtbasetbl db 0,1,4,9,16,25,36,49,64,81,100,121,144,169,196,225

           public sqrt

           align 16

sqrt:
           pushad
           mov ebp,eax
           bsr ebx,eax
           jnz short sqrtf0
           xor ebx,ebx
sqrtf0:
           shr ebx,3
           lea eax,[ebx*8]
           mov cl,32
           sub cl,al
           rol ebp,cl
           mov eax,ebp
           movzx eax,al
           mov edi,offset sqrtbasetbl
           mov ecx,10h
sqrtl0:
           scasb
           je short sqrtl0d
           jb short sqrtl0d2
           loop sqrtl0
           inc edi
sqrtl0d2:
           dec edi
           inc cl
sqrtl0d:
           movzx edx,byte ptr [edi-1]
           dec cl
           xor cl,0fh
           mov edi,ecx
           mov ecx,ebx
           jecxz short sqrtdone
           sub eax,edx
sqrtml:
           shld eax,ebp,8
           rol ebp,8
           mov ebx,edi
           shl ebx,5
           xor edx,edx
           mov esi,eax
           div ebx
           rol edi,4
           add edi,eax
           add ebx,eax
sqrtf2:
           imul eax,ebx
           mov edx,eax
           mov eax,esi
           sub eax,edx
           jc short sqrtf1
           loop sqrtml
sqrtdone:
           mov [esp+28],edi
           popad
           ret
sqrtf1:
           dec ebx
           dec edi
           movzx eax,bl
           and al,1fh
           jmp sqrtf2

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Set_precal147: Set precal7 for plane transformation
;
; In:
;   ECX - y location to pre-calculate plane to
; Out:
;   precal1
;   precal4
;   precal7
;
; Notes:
; Plane is ecx and allows above formulas to determine where a
; point/object would be along that plane if z is not negative
;
; Good for runway translations or super huge background polygons - not used
; by regular 3d.asm routines
;
; How to use: lets say you've got a billion background objects that are
; on the ground (or all on the same y plane).  you call set_precal147
; with that y plane location and use frotate instead of erotate to
; determine where points rotated along that plane will end up.  this
; speeds the routine up by 33% by cutting out 3 imuls.
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

            align 16

precal1     dd 0
precal4     dd 0
precal7     dd 0

set_precal147:
            if usez eq yes
            mov eax,ecx
            sub eax,eyey
            imul ematrix+4
            shrd eax,edx,14
            mov precal1,eax
            endif

            mov eax,ecx
            sub eax,eyey
            imul ematrix+16
            shrd eax,edx,14
            mov precal4,eax

set_precal7:
            sub ecx,eyey
            mov eax,ecx
            imul ematrix+28
            shrd eax,edx,14
            mov precal7,eax
            ret

            align 16

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Frotate: fast object/point rotation along pre-calculated y plane
;
; In:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;    ematrix - 32 bit rotation matrix - set up by "setsincose" routine
;    precal1 - set up by set_precal147
;    precal4 - set up by set_precal147
;    precal7 - set up by set_precal147
;
; Out:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;
; Notes:
;    see above
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

frotate:
           mov eax,ematrix+8
           imul ebp
           shrd eax,edx,14
           mov edi,eax
           if usez eq yes
           add edi,precal1
           endif
           mov eax,ematrix+0
           imul ebx
           shrd eax,edx,14
           add edi,eax   ; di = new x

           mov eax,ematrix+20
           imul ebp
           shrd eax,edx,14
           mov esi,eax
           add esi,precal4
           mov eax,ematrix+12
           imul ebx
           shrd eax,edx,14
           add esi,eax   ; si = new y

           mov eax,ematrix+32
           imul ebp
           shrd eax,edx,14
           mov ebp,eax
           add ebp,precal7
           mov eax,ematrix+24
           imul ebx
           shrd eax,edx,14
           add ebp,eax   ; bp = new z

           mov ecx,esi
           mov ebx,edi

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Fzrotate: fast single variable solve for object/point rotation along
;           pre-calculated y plane
;
; In:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;    ematrix - 32 bit rotation matrix - set up by "setsincose" routine
;    precal1 - set up by set_precal147
;    precal4 - set up by set_precal147
;    precal7 - set up by set_precal147
;
; Out:
;    ESI - z point
;
; Notes:
; Fast solve for single matrix variable similar to erotate but uses frotate
; plane matrix with precal147
;
; remember , matrix offsets are:
;
;  0 1 2     multiply those by 4 for the doublewords
;  3 4 5
;  6 7 8
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           align 16

fzsolve:
           mov eax,ematrix+32      ; solve z = bx(6)+cx(7)+bp(8)
           imul ebp
           shrd eax,edx,14
           mov esi,eax
           add esi,precal7
           mov eax,ematrix+24
           imul ebx
           shrd eax,edx,14
           add esi,eax   ; si = new z
           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Fxrotate: fast single variable solve for object/point rotation along
;           pre-calculated y plane
;
; In:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;    ematrix - 32 bit rotation matrix - set up by "setsincose" routine
;    precal1 - set up by set_precal147
;    precal4 - set up by set_precal147
;    precal7 - set up by set_precal147
;
; Out:
;    EDI - x point
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           align 16
fxsolve:
           mov eax,ematrix+8        ; solve x = bx(0)+cx(1)+bp(2)
           imul ebp
           shrd eax,edx,14
           mov edi,eax
           if usez eq yes
           add edi,precal1
           endif
           mov eax,ematrix+0
           imul ebx
           shrd eax,edx,14
           add edi,eax   ; di = new x
           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Fyrotate: fast single variable solve for object/point rotation along
;           pre-calculated y plane
;
; In:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;    ematrix - 32 bit rotation matrix - set up by "setsincose" routine
;    precal1 - set up by set_precal147
;    precal4 - set up by set_precal147
;    precal7 - set up by set_precal147
;    ESI - z point
;    EDI - x point
; Out:
;    EBX - x point
;    ECX - y point
;    EBP - z point
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           align 16
fysolve:
           mov eax,ematrix+20       ; solve y = bx(3)+cx(4)+bp(5)
           imul ebp
           shrd eax,edx,14
           mov ecx,eax
           add ecx,precal4
           mov eax,ematrix+12
           imul ebx
           shrd eax,edx,14
           add ecx,eax   ; cx = new y

           mov ebx,edi
           mov ebp,esi

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Lambert: generate lambert shading 1x3 matrix, completed in 6 imuls
;
; In:
;    ESI - Object # to get angles from
;    vxs[esi*2] - object x angle (0-65536)
;    vys[esi*2] - object y angle (0-65536)
;    vzs[esi*2] - object z angle (0-65536)
;    y_angle_of_sun - (0-65536)
;
; Out:
;    lmatrix  - shading matrix
;    ESI - ?
;
; Notes:
;
;z= ( sz ( cx + ( sx * cy )) + cz * sy ) * 45degrees  [x]
;   ( cz ( cx + ( sx * cy )) - sz * sy ) * 45degrees  [y]
;   ( cx * cy - sx ) * 45 degrees                     [z]
;
;note cos45=sin45=2d41h, but we will use 2d00h (99.2% accurate)
; you can change the y angle of the sun/light but not the x angle.
; changing the x angle would require a new formula.
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

lambert:
           mov ax,vxs[esi*2]
           neg ax
           push eax
           call cosign
           mov vcosx,eax
           pop eax
           call sign
           mov vsinx,eax
           mov ebp,eax            ; ebp = sx

           mov ax,vzs[esi*2]
           neg ax
           push eax
           call cosign
           mov vcosz,eax
           mov edi,eax            ; edi = cz
           pop eax
           call sign
           mov vsinz,eax
           mov edx,eax            ; edx = sz

           mov ax,vys[esi*2]
           neg ax
           add eax,y_angle_of_sun ; 2000h = 45 degrees y angle for light source
           push eax
           call cosign
           mov vcosy,eax
           mov esi,eax            ; esi = cy
           pop eax
           call sign
           mov vsiny,eax          ; eax = sy

           mov ebx,edx            ; ebx = sz

           mov ecx,eax            ; ecx = sy

           mov eax,ebp            ; get sx

           imul esi               ; eax = sx * cy
           shrd eax,edx,14
           sub eax,vcosx          ; eax = cx + ( sx * cy)

           push eax

           imul ebx               ; cx + ( sx * cy) * sz
           shrd eax,edx,14
           mov lmatrix+0,eax

           pop eax

           imul edi               ; di = cz
           shrd eax,edx,14
           mov lmatrix+4,eax      ; cx + ( sx * cy) * cz

           mov eax,ebx
           imul ecx               ; - sz * sy
           shrd eax,edx,14
           sub lmatrix+4,eax

           mov eax,edi
           imul ecx               ; cz * sy
           shrd eax,edx,14
           add lmatrix+0,eax

           mov eax,vcosx          ; (cx * cy - sx) * 45deg
           imul esi
           shrd eax,edx,14
           mov ebx,eax
           add ebx,ebp

           cmul eax,ebx,2d00h     ; * 45degrees
           shrd eax,edx,14
           movsx eax,ax
           mov lmatrix+8,eax

           mov ebx,lmatrix+4
           cmul eax,ebx,2d00h
           shrd eax,edx,14
           mov lmatrix+4,eax

           mov ebx,lmatrix+0
           cmul eax,ebx,2d00h
           shrd eax,edx,14
           mov lmatrix+0,eax

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Pre_cal_lambert: Pre-calculate all neccessary stuff for object DI
;
; In:
;    EDI - Object # to pre-calculate normals for
;    objbase[esi*4] -> offset of object data
; Out:
;    ESI -> minimum address of object data
;    EDI -> maximum address of object data
;    EBP -> points to header for object
;
; Notes:
;
; Precalculate surface normals for object di.  This  is  so  you  don't
; have to type them in when designing new objects. Imagine, 400 points,
; with 350 surfaces, calculating them all manually?  This routine  also
; figures out the iteration skip offset (if you have surfaces dependant
; on other surfaces) and also sets bit 1 if it is a line (two  points),
; and sets bit 4 if 1 point.This routine also sets the number of points
; to skip if an iteration is found. It  counts  the  number  of  points
; within iterations (even iterations within iterations)  and  sets  the
; skip value so any iterations skipped will have a pre-calculated point
; offset.  Did that make sense?
;
; Things done here:
;
; set point command if only 1 connection
; set line  command if only 2 connections
; set normal bit in commands if shading used in texture
; calculate and set shading normals
; calculate offsets for iteration jumps (in case surface not visible)
; calculate number of points to skip for iterations (in case surface not visible)
; set offset flag if iteration uses a point offset (4'th future use word)
; calculate and set auto-intensity of color if auto_s bit set
;
; Most of the above is done so the user (you) wont have to calculate this stuff
; yourself - makes object modification much easier.
;
; If you find the routine to be sloppy remember it is only used
; for object initialization.
;
; This routine  will probably crash if your object is not set up correctly
; The entire 3dvect source will crash if this routine  isn't run  and  the
; chances of you knowing how to do all that  this  routine  does  manually
; are pretty slim since most of you out there are doughheads.
;
; The minimum and maximum addresses are returned in ESI and EDI so you can
; output the resulting object to a binary file. EBP points to the starting
; header of the object (usually, but not neccessaraly, the minimum address)
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

lx1  dd 0
ly1  dd 0
lz1  dd 0

lx2  dd 0
ly2  dd 0
lz2  dd 0

lx3  dd 0
ly3  dd 0
lz3  dd 0

finx dd 0
finy dd 0
finz dd 0

temp1 dw 0
temp2 dw 0 ; number of points
temp3 dw 0 ; number of sides
temp4 dw 0
temp5 dd 0 ; minimum address
temp6 dd 0 ; maximum address

pre_cal_lambert:
           movzx edi,di           ; in case user is lazy
           mov esi,objbase[edi*4]
           mov temp5,esi
           mov temp6,esi
           push esi               ; save for exit
more_reses:
           push esi               ; save header offset
           add esi,4
           lodsd
           add esi,eax            ; handle first resolution

           lodsw
           mov temp2,ax
           lodsw
           mov temp3,ax

           mov xad,0
           mov yad,0
           mov zad,0
           mov temp1,-1
           mov ax,[esi+14]
           mov temp4,ax

           mov eax,[esi+8*2]
           or eax,[esi+8*2+4]
           jnz lam_hhgg
           mov dword ptr [esi+8*2],maxz    ; if no max/min found, force one
           mov dword ptr [esi+8*2+4],minz
lam_hhgg:
           mov eax,[esi+12*2]
           jnz lam_hhgc
           mov dword ptr [esi+12*2],tolerance
lam_hhgc:
           add esi,25*2           ; skip future use bytes

           mov edi,4              ; edi=2 to skip center of gravity
           mov xp,0
           mov yp,0
           mov zp,0
           cmp temp2,0
           je no_points_2
lam_ap12:
           mov bx,w [esi]         ; load all the points into array
           mov cx,w [esi+2]
           mov bp,w [esi+4]
           add bx,w xad
           add cx,w yad
           add bp,w zad
           movsx ebx,bx
           movsx ecx,cx
           movsx ebp,bp
           mov xp[edi],ebx
           mov yp[edi],ecx
           mov zp[edi],ebp
           add esi,6
           add edi,4
           dec temp2
           jne s lam_ap12         ; esi = address of sides now...
no_points_2:
           mov pointindex,edi

lam_loadsides:
           call checkesi          ; check minimum and maximum addresses

           mov edi,esi            ; save in case of line adjust
           mov ax,[esi]           ; get command

           mov bx,ax              ; save command
           test ax,special        ; test if special command
           jz s lam_notmap        ; no, skip through loop
           mov di,ax
           and edi,special-1

           cmp ax,gosub           ; check for jump commands, yeah,yeah, I should have made a table, who cares...
           jne ffgg1
           add esi,2
           lodsw
           push esi
           movsx eax,ax
           sub esi,2
           add esi,eax
           jmp lam_next
ffgg1:
           cmp ax,return
           jne ffgg2
           pop esi
           jmp lam_next
ffgg2:
           cmp ax,goto_offset
           jne ffgg3
           add esi,2
           lodsw
           movsx eax,ax
           sub esi,2
           add esi,eax
           jmp lam_next
ffgg3:
           mov di,number_ofb[edi*2] ; yes, skip special command length
           add esi,edi
           mov bx,0
           cmp ax,sub_object
           je lam_do_it
           cmp ax,static_sub_object
           je lam_do_it
           jmp lam_next           ; go to next side
lam_notmap:
           call pcl_testashade

           mov ax,[esi+2]         ; get texture for both sides
           or  ax,[esi+4]

           test ax,shade          ; test shading bit
           jnz lam_calcit         ; yes, calculate shading normal

           push edi               ; save command location

           add esi,4+4+2          ; skip 2 colour & 2 texture words & command
           mov edi,esi

           lodsw                  ; get first point indexer
           add ax,temp4
           stosw

           mov cx,ax
           mov dx,0
lam_ldlp:
           lodsw                  ; count number of connection points
           add ax,temp4
           stosw

           inc dx
           cmp ax,cx
           jne lam_ldlp

           call checkesi          ; check minimum and maximum addresses

           pop edi                ; pop command location

           cmp dx,1               ; only 1 point?, set point command
           jne lam_test_line
           or  w [edi+0],both
           or  w [edi+2],point
           or  w [edi+4],point
           jmp lam_test_iteration

lam_test_line:
           cmp dx,2               ; only 2 points?, set line command
           jne lam_test_iteration
           or  w [edi+0],both
           or  w [edi+2],line
           or  w [edi+4],line

lam_test_iteration:
           mov ax,0
           test bx,iterate        ; test if iteration command used
           jnz lam_do_it          ; yes,solve internal iteration
lam_next:
           dec temp1
           jnz lam_nopop
           pop ax
           mov temp4,ax
           pop ax
           mov temp1,ax
           jmp lam_next
lam_nopop:
           dec temp3
           jnz lam_loadsides

           pop esi
           lodsd
           add esi,4
           cmp eax,-1             ; last resolution?
           jne more_reses

           call checkesi          ; check minimum and maximum addresses
           pop ebp                ; pop header offset
           mov esi,temp5          ; load up start and end offsets of object
           mov edi,temp6

           ret

lam_calcit:
           push esi               ; save command location
           add esi,4+4+2          ; skip colour and 2 future use words

           lodsw                  ; first point
           add ax,temp4
           mov [esi-2],ax
           push ax
           movzx edi,ax
           shl edi,2

           mov ebx,[xp+edi]
           mov ecx,[yp+edi]
           mov ebp,[zp+edi]

           mov lx1,ebx
           mov ly1,ecx
           mov lz1,ebp

           lodsw                  ; second point
           add ax,temp4
           mov [esi-2],ax
           movzx edi,ax
           shl edi,2

           mov ebx,[xp+edi]
           mov ecx,[yp+edi]
           mov ebp,[zp+edi]

           mov lx2,ebx
           mov ly2,ecx
           mov lz2,ebp

           lodsw                  ; third point
           add ax,temp4
           mov [esi-2],ax
           movzx edi,ax
           shl edi,2

           mov ebx,[xp+edi]
           mov ecx,[yp+edi]
           mov ebp,[zp+edi]

           mov lx3,ebx
           mov ly3,ecx
           mov lz3,ebp

           push esi

           call calc_normal

           pop esi

           pop dx                 ; now find shading normal storage, pop first connector

lam_ldl2:
           lodsw
           add ax,temp4
           mov [esi-2],ax
           cmp ax,dx
           jne lam_ldl2

           mov edi,esi

           mov ax,bx
           stosw
           mov ax,cx
           stosw
           mov ax,bp
           stosw

           add esi,6

           pop edi         ; get original command location back
           or w [edi],normal
           mov bx,[edi]
           jmp lam_test_iteration

lam_surfc_cnt dw 0

; this finds the total number of points to skip if an iteration fails, dx = #
; remember, this is a pre-calculation routine so it doesn't need to be fast.

lam_do_it:
           mov dx,0        ; clear total number of points to skip

           mov ax,[esi+10]        ; test if there is a center of rotation point
           or ax,[esi+12]
           or ax,[esi+14]
           jz done_alter2
           or w [esi+8],centroid  ; set flag if offset found (center of gravity)

           movsx ebx,w [esi+10]
           movsx ecx,w [esi+12]
           movsx ebp,w [esi+14]
           add xad,ebx
           add yad,ecx
           add zad,ebp
           mov ebx,xad
           mov ecx,yad
           mov ebp,zad
           mov edi,pointindex
           mov xp[edi],ebx
           mov yp[edi],ecx
           mov zp[edi],ebp
           add pointindex,2
           add dx,1        ; centroid is an extra point to skip
done_alter2:
           mov ax,temp1
           push ax
           mov ax,temp4
           push ax

           mov ax,[esi+18]
           add temp4,ax
           mov w [esi+18],0

           mov lam_surfc_cnt,0
           push esi        ; this is our return address (continue from here+4)

           lodsw           ; get number of points.
           add dx,ax       ; save as TOTAL number of points to skip
           mov temp2,ax

           lodsw           ; get number of surfaces
           mov temp1,ax
           inc temp1
           add lam_surfc_cnt,ax ; count until this is zero

           mov eax,[esi+8*2]
           or eax,[esi+10*2]
           jnz lam_hhgg2
           mov dword ptr [esi+8*2],maxz  ; if no max/min found, force one
           mov dword ptr [esi+8*2+4],minz
lam_hhgg2:
           mov eax,[esi+12*2]
           jnz lam_hhgl
           mov dword ptr [esi+12*2],tolerance
lam_hhgl:
           add esi,25*2

           mov edi,pointindex
           cmp temp2,0
           je lam_test_check      ; only sides added, no additional points
lam_ap13:
           mov bx,w [esi]         ; load all the points into array
           mov cx,w [esi+2]       ; for calculation of gourad shadings
           mov bp,w [esi+4]
           add bx,w xad
           add cx,w yad
           add bp,w zad
           movsx ebx,bx
           movsx ecx,cx
           movsx ebp,bp
           mov xp[edi],ebx
           mov yp[edi],ecx
           mov zp[edi],ebp
           add esi,6
           add edi,4
           dec temp2
           jne s lam_ap13         ; esi = address of sides now...

           mov pointindex,edi

lam_test_check:
           cmp lam_surfc_cnt,0    ; test if user just wants to add points
           je lam_no_surfs        ; i dont know why anyone would want to do this?

lam_test_until_target:
           lodsw                  ; get command
           mov bx,ax
           test ax,special        ; test if special command
           jz s lam_notmap_it     ; no, skip through loop

           mov di,ax
           and edi,special-1

           cmp ax,gosub           ; check for jump commands, yeah,yeah, I should have made a table, who cares...
           jne ffgg1x
           add esi,2
           lodsw
           push esi
           movsx eax,ax
           sub esi,2
           add esi,eax
           jmp lam_next_it
ffgg1x:
           cmp ax,return
           jne ffgg2x
           pop esi
           jmp lam_next_it
ffgg2x:
           cmp ax,goto_offset
           jne ffgg3x
           add esi,2
           lodsw
           movsx eax,ax
           sub esi,2
           add esi,eax
           jmp lam_next_it
ffgg3x:
           movzx edi,number_ofb[edi*2] ; yes, skip special command length
           add esi,edi
           sub esi,2
           cmp ax,sub_object
           je lam_re_lam
           cmp ax,static_sub_object
           je lam_re_lam
           jmp s lam_nog

lam_notmap_it:
           lodsw
           mov bp,ax
           lodsw
           or bp,ax               ; find if shading bit used, add esi,6 if so
           add esi,4              ; skip 2 colour words

           lodsw                  ; get first point indexer
           mov cx,ax
lam_ldl3:
           lodsw
           cmp ax,cx
           jne lam_ldl3

           test bp,shade          ; test if gouraud normal present
           jz lam_nog
           add esi,6              ; skip it if present
lam_nog:
           test bx,iterate        ; test if iteration command used
           jnz lam_re_lam         ; solve internal iteration again...
lam_next_it:
           dec lam_surfc_cnt
           jnz lam_test_until_target
lam_no_surfs:
           mov edi,esi            ; save current location
           pop esi                ; return original start location
           sub edi,esi            ; get difference between them
           sub di,8               ; offset for loadsides routine - constant

           lodsw                  ; get number of points
           mov bx,ax

           lodsw
           add temp3,ax

           mov cx,dx
           mov ax,6
           imul bx
           movzx ebx,ax

           mov ax,di
           mov edi,esi
           stosw                  ; save offset
           mov ax,cx
           stosw                  ; save number of points found in iterations

           mov eax,[esi+8*2]
           or eax,[esi+10*2]
           jnz lam_hhgg3
           mov dword ptr [esi+8*2],maxz  ; if no max/min found, force one
           mov dword ptr [esi+8*2+4],minz
lam_hhgg3:
           mov eax,[esi+12*2]
           jnz lam_hhgq
           mov dword ptr [esi+12*2],tolerance
lam_hhgq:
           add esi,25*2           ; adjust for next load
           add esi,ebx

           jmp lam_next

lam_re_lam:
           lodsw           ; get number of points for recursed iteration
           add dx,ax       ; save as TOTAL number of points to skip
           mov cx,ax

           lodsw           ; get number of surfaces
           add lam_surfc_cnt,ax ; count until this is zero

           mov eax,[esi+8*2]
           or eax,[esi+10*2]
           jnz lam_hhgg4
           mov dword ptr [esi+8*2],maxz  ; if no max/min found, force one
           mov dword ptr [esi+8*2+4],minz
lam_hhgg4:
           mov eax,[esi+12*2]
           jnz lam_hhgv
           mov dword ptr [esi+12*2],tolerance
lam_hhgv:
           add esi,25*2

           mov eax,6
           imul cx
           add esi,eax

           jmp lam_next_it

checkesi:
           cmp esi,temp5
           jae noesi1
           mov temp5,esi
noesi1:
           cmp esi,temp6
           jbe noesi2
           mov temp6,esi
noesi2:
           ret

pcl_testashade:
           mov ax,[esi+2]
           or  ax,[esi+4]
           test ax,auto_s
           jnz pcl_test2
           ret
pcl_test2:
           pushad
           push esi
           push esi
           mov vxs,0
           mov vys,0
           mov vzs,0
           mov esi,0
           call lambert
           pop esi

           movzx ebx,w [esi+10]
           movzx ecx,w [esi+12]
           movzx edx,w [esi+14]

           mov eax,xp[ebx]
           mov lx1,eax
           mov eax,yp[ebx]
           mov ly1,eax
           mov eax,zp[ebx]
           mov lz1,eax
           mov eax,xp[ecx]
           mov lx2,eax
           mov eax,yp[ecx]
           mov ly2,eax
           mov eax,zp[ecx]
           mov lz2,eax
           mov eax,xp[edx]
           mov lx3,eax
           mov eax,yp[edx]
           mov ly3,eax
           mov eax,zp[edx]
           mov lz3,eax
           call calc_normal

           call lrotate
           pop esi
           push edi

           test w [esi+2],auto_s
           jz pcl_test4
           test w [esi+2],inverse
           jz pcl_test3
           neg edi
pcl_test3:
           add edi,256
           shr edi,1              ; result -256 to +256, turn into 0-256
           mov al,b shading_tables[edi] ; now into 0-15
           xor ah,ah
           add w [esi+6],ax
           mov ax,[esi+2]
           and ax,65535-auto_s
           mov [esi+2],ax
pcl_test4:
           pop edi

           test w [esi+4],auto_s
           jz pcl_test6
           test w [esi+4],inverse
           jz pcl_test5
           neg edi
pcl_test5:
           add edi,256
           shr edi,1              ; result -256 to +256, turn into 0-256
           mov al,b shading_tables[edi] ; now into 0-15
           xor ah,ah
           add w [esi+8],ax
           mov ax,[esi+4]
           and ax,65535-auto_s
           mov [esi+4],ax
pcl_test6:
           popad
           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Calc_normal: calculate surface normal
;
; In:
;    LX1 - x of point 1 on triangle
;    LY1 - y of point 1 on triangle
;    LZ1 - z of point 1 on triangle
;    LX2 - x of point 2 on triangle
;    LY2 - y of point 2 on triangle
;    LZ2 - z of point 2 on triangle
;    LX3 - x of point 3 on triangle
;    LY3 - y of point 3 on triangle
;    LZ3 - z of point 3 on triangle
;
; Out:
;    EBX = finx = x of surface normal of triangle
;    ECX = finy = y of surface normal of triangle
;    EBP = finz = z of surface normal of triangle
;
; Notes:
; x2 = x2 - x1
; y2 = y2 - y1
; z2 = z2 - z1
;
; x3 = x3 - x1
; y3 = y3 - y1
; z3 = z3 - z1
;
; x = y2 * z3 - z2 * y3
; y = z2 * x3 - x2 * z3
; z = x2 * y3 - y2 * x3
;
; a = SQR(x ^ 2 + y ^ 2 + z ^ 2)
;
; x = INT(x / a * 256 + .5)
; y = INT(y / a * 256 + .5)
; z = INT(z / a * 256 + .5)
;
; This worked for me on the first try!
;
; If you wanted to get the equation of a plane, you could do this after:
;  d = - x * x1 - y * y1 - z * z1
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

           nshl = 8

calc_normal:
           mov ebx,lx1
           mov ecx,ly1
           mov ebp,lz1

           sub lx2,ebx
           sub ly2,ecx
           sub lz2,ebp

           sub lx3,ebx
           sub ly3,ecx
           sub lz3,ebp

           mov eax,ly2
           mov ebx,lz3
           imul ebx
           mov ecx,eax

           mov eax,lz2
           mov ebx,ly3
           imul ebx
           sub ecx,eax

           mov finx,ecx ; save x of normal

           mov eax,lz2
           mov ebx,lx3
           imul ebx
           mov ecx,eax

           mov eax,lx2
           mov ebx,lz3
           imul ebx
           sub ecx,eax

           mov finy,ecx ; save y of normal

           mov eax,lx2
           mov ebx,ly3
           imul ebx
           mov ecx,eax

           mov eax,ly2
           mov ebx,lx3
           imul ebx
           sub ecx,eax

           mov finz,ecx ; save z of normal

calc_testloop:
           cmp finx,32768 ; make sure (normal^2)*2 is < 2^32
           jge calc_shrit
           cmp finy,32768
           jge calc_shrit
           cmp finz,32768
           jge calc_shrit

           cmp finx,-32768
           jle calc_shrit
           cmp finy,-32768
           jle calc_shrit
           cmp finz,-32768
           jg  ok_2_bite_dust

calc_shrit:
           shr finx,1   ; calculations will be too large if squared, div by 2
           test finx,40000000h
           jz no_neg_calc1
           or   finx,80000000h
no_neg_calc1:
           shr finy,1
           test finy,40000000h
           jz no_neg_calc2
           or   finy,80000000h
no_neg_calc2:
           shr finz,1
           test finz,40000000h
           jz no_neg_calc3
           or   finz,80000000h
no_neg_calc3:
           jmp calc_testloop

ok_2_bite_dust:
           mov eax,finx ; x^2
           mov edi,eax  ; objects
           imul edi
           mov edi,eax

           mov eax,finy ; y^2
           mov esi,eax
           imul esi
           mov esi,eax

           mov eax,finz ; z^2
           mov ebp,eax
           imul ebp

           add eax,esi
           add eax,edi

           call sqrt    ; get square root of number

           mov ecx,eax
           cmp ecx,0
           je lam_abort ; should never happen!

           mov eax,finx
           cdq
           shl eax,nshl ; set unit vector to 2^nshl (256)
           idiv ecx
           mov finx,eax

           mov eax,finy
           cdq
           shl eax,nshl
           idiv ecx
           mov finy,eax

           mov eax,finz
           cdq
           shl eax,nshl
           idiv ecx
           mov finz,eax

           mov ebx,finx
           mov ecx,finy
           mov ebp,finz

lam_abort:
           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
; Calc_D: Calculate D portion of equation of a plane
; In:
;    EBX = x of surface normal of triangle
;    ECX = y of surface normal of triangle
;    EBP = z of surface normal of triangle
;    LX1 - x of point on triangle (any point)
;    LY1 - y of point on triangle
;    LZ1 - z of point on triangle
; Out:
;    EAX = D     (Ax+By+Cz=D)
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
calc_d:
           mov eax,lx1
           imul ebx
           mov esi,eax

           mov eax,ly1
           imul ecx
           add esi,eax

           mov eax,lz1
           imul ebp
           add eax,esi

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Set_up_all_lambert: set up all lambert normals from object si to object di
;
; In:
;    ESI - object # to start at
;    EDI - object # to end at
;    objbase[esi*4 - edi*4] -> offsets to object data
; Out:
;   null
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

set_up_all_lambert:
           movzx edi,di ; in case user is lazy
           movzx esi,si

           xchg edi,esi ; so user doesn't get confuzed
set_lop:
           push esi edi
           call pre_cal_lambert
           pop edi esi
           inc edi
           cmp edi,esi
           jna  set_lop

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Lrotate: rotate surface normal through lambert matrix
;
; In:
;    BX - x of surface normal
;    CX - y of surface normal
;    BP - z of surface normal
;    lmatrix - 16 bit, 1x3 lambert shading matrix - set up by "lambert" routine
; Out:
;    BX - x of surface normal (untouched)
;    CX - y of surface normal (untouched)
;    BP - z of surface normal (untouched)
;   EDI - colour intensity for surface (-255 to +255)
;
; Notes:
;   Your mother is a hamster.
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

lrotate:
           mov eax,lmatrix+8    ; solve edi = bx(0)+cx(4)+bp(8)
           imul bp
           shrd ax,dx,14
           mov edi,eax
           mov eax,lmatrix+4
           imul cx
           shrd ax,dx,14
           add edi,eax
           mov eax,lmatrix+0
           imul bx
           shrd ax,dx,14
           add edi,eax          ; di = new colour -256 to 255 (not edi, di!)
           movsx edi,di

           add edi,256          ; make sure result is within range -256 to 256
           and edi,511          ; sometimes is messes up...
           sub edi,256

           ret

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Shading tables used for fake cosin colour intensity - 256 bytes
; Default is 16 colours per lambert calculation.  But you could have 32, 48
; or whatever you want, even an odd number like 53.  Use the SHADING.BAS
; program to make the table to your custom size.
;
; Shading_bits is the variable for use with the texture command "LAST".  This
; variable tells the routine what bits to  pluck  off  when  determining  the
; shading intensity. Obviosly this cant  be  used  if  your  palette  shading
; length is not a function of 2's complement.
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

  if shading_colours eq 16

; 16 colour shading table

  shading_bits equ 0fh

  align 16
shading_tables:
  db 0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1
  db 1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2
  db 2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3
  db 3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4
  db 4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5
  db 5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6
  db 6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7
  db 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7
  db 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8
  db 8,8,8,9,9,9,9,9,9,9,9,9,9,9,9,9
  db 9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10
  db 10,10,10,10,10,10,11,11,11,11,11,11,11,11,11,11
  db 11,11,11,11,11,11,11,11,12,12,12,12,12,12,12,12
  db 12,12,12,12,12,12,12,12,12,13,13,13,13,13,13,13
  db 13,13,13,13,13,13,13,13,13,14,14,14,14,14,14,14
  db 14,14,14,14,14,14,14,14,15,15,15,15,15,15,15,15

  elseif shading_colours eq 32

; 32 colour shading table

  shading_bits equ 1fh

   align 16
shading_tables:
   db 0,0,0,0,0,0,0,1,1,1,1,1,1,1,2,2
   db 2,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4
   db 4,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6
   db 6,6,6,6,7,7,7,7,7,7,7,7,8,8,8,8
   db 8,8,8,8,9,9,9,9,9,9,9,9,10,10,10,10
   db 10,10,10,10,10,11,11,11,11,11,11,11,11,12,12,12
   db 12,12,12,12,12,12,13,13,13,13,13,13,13,13,14,14
   db 14,14,14,14,14,14,14,15,15,15,15,15,15,15,15,15
   db 16,16,16,16,16,16,16,16,16,17,17,17,17,17,17,17
   db 17,18,18,18,18,18,18,18,18,18,19,19,19,19,19,19
   db 19,19,20,20,20,20,20,20,20,20,20,21,21,21,21,21
   db 21,21,21,22,22,22,22,22,22,22,22,22,23,23,23,23
   db 23,23,23,23,24,24,24,24,24,24,24,24,25,25,25,25
   db 25,25,25,25,26,26,26,26,26,26,26,26,27,27,27,27
   db 27,27,27,28,28,28,28,28,28,28,29,29,29,29,29,29
   db 29,29,30,30,30,30,30,30,30,31,31,31,31,31,31,31
   endif

;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░
;
; Sqrax2bx2: Fast pathagorean solve, sqr(ax^2+bx^2)
;
; In:
;   EAX - a number
;   EBX - another number
; Out:
;   EAX - sqr(EAX^2+EBX^2)
;
; Notes:
;
; Works well for numbers eax<90, ebx<128.  Uses huge table.
; Limit of solve is for 32 bit values. Only works for positive values of ax&bx
; Numbers outside the limit of the tabel could be shr'd and re-tested rather
; than using the mathematical solution. There should be a .BAS Qbasic file
; that generates the table provided with the 3Dvect source.
;
;░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░░

pathagorean:
;          include pathagor.inc ; dont include this routine if you don't need it
                                ; 'cause it's a real big table.
sqrax2bx2:
           cmp eax,ebx          ; set eax = smallest
           ja s pa_smal
           xchg eax,ebx
pa_smal:
           cmp eax,90           ; check if parameters of triangle are within
           jae s pa_imul        ; table or are too big for fast load.
           cmp ebx,128
           jae s pa_imul

           shl ebx,7            ; *128
           mov al,b pathagorean[eax+ebx]
           ret
pa_imul:
           imul eax,eax         ; variables too large, do mathematically
           imul ebx,ebx
           add eax,ebx
           jmp sqrt